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Abstract 

This  report  summarizes  the  research  objectives  achieved  in  this  project  during  the  period  03- 
01-2009  to  12-31-2011.  Computational  techniques  have  been  developed  that  enable  electronic 
structure  calculations  at  macroscopic  scales  by  seamless  bridging  of  the  quantum  length-scale 
with  continuum.  The  various  components  of  the  developed  multi-scale  scheme  include  (i) 
real-space  formulation  of  density-functional  theory  (DFT);  (ii)  finite-element  discretization; 
(iii)  an  adaptive  coarse-graining  through  quasi-continuum  reduction.  The  formulation  has 
been  developed  for  orbital-free  DFT  with  non-local  kinetic  energy  functionals.  Studies  on 
the  energetics  of  defects  in  materials  using  the  developed  techniques  have  provided  many 
new  insights  into  the  behavior  of  defects,  and  the  complex  nature  of  the  interacting  length- 
scales  that  influence  their  behavior.  An  extension  of  this  formulation  to  Kohn-Sham  DFT 
was  attempted,  and  as  part  of  this  effort  an  efficient  real-space  formulation  of  Kohn-Sham 
DFT  was  developed.  Studies  show  that  the  developed  formulation  compares  favorably  with 
existing  conventional  DFT  implementations,  enables  consideration  of  complex  geometries 
and  boundary  conditions,  and  exhibits  good  scalability  on  parallel  computing  architectures. 


Summary  of  Research  Objectives  Achieved 

In  this  project  we  have  conducted  a  program  of  research  to  develop  computational  tech¬ 
niques  that  seamlessly  bridge  the  quantum  length-scale  with  the  continuum  and  enable 
an  accurate  electronic  structure  study  of  defective  crystals.  The  techniques  have  been  de¬ 
veloped  for  orbital- free  DFT  (Wang  &  Teter,  1992;  Wang  et  al.,  1998,  1999),  which  is  an 
approximation  to  the  Kohn-Sham  density  functional  theory  where  the  kinetic  energy  func¬ 
tionals  are  modeled  and  is  applicable  to  materials  systems  whose  electronic  structure  is  close 
to  a  free  electron  gas.  Further,  in  order  to  extend  these  techniques  to  Kohn-Sham  den¬ 
sity  functional  theory,  we  have  developed  a  real-space  formulation  of  Kohn-Sham  density 
functional  theory  (Suryanarayana  et  ah,  2010;  Motamarri  et  ah,  2012),  which  is  the  basis 
for  the  coarse-graining  techniques  based  on  quasi-continuum  reduction  (Gavini  et  ah,  2007; 
Radhakrishnan  &  Gavini,  2010).  A  comprehensive  summary  of  the  important  research  ob¬ 
jectives  achieved  in  this  project  is  presented  below  with  references  to  the  published  articles. 
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1  Development  of  quasi-continuum  reduction  of  orbital- 
free  DFT  with  non-local  kinetic  energy  functionals 
and  Mathematical  Analysis: 

We  have  accomplished  the  quasi- continuum  reduction  of  orbital-free  DFT  with  non-local 
kinetic  energy  functionals  (Radhakrishnan  &  Gavini,  2010),  which  has  enabled  electronic 
structure  calculations  using  orbital-free  DFT  without  cell-size  restrictions  and  has  enabled 
an  accurate  electronic  structure  study  of  defects  within  the  framework  of  orbital-free  DFT. 
One  of  the  significant  challenges  in  the  development  of  the  formulation  was  the  incorpora¬ 
tion  of  non-local  kinetic  energy  functionals,  density-independent  (DI)  and  density-dependent 
(DD)  functionals,  into  the  quasi-continuum  reduction  which  is  based  on  a  local  variational 
formulation.  These  non-local  kernel  energies  constitute  an  important  part  of  the  orbital-free 
kinetic  energy  functional,  which  are  necessary  to  satisfy  the  Lindhard  response  function — a 
known  limiting  behavior  of  uniform  electron  gas.  The  Wang-Govind-Carter  (WGC)  kinetic 
energy  functional  (Wang  et  ah,  1999),  which  falls  into  the  family  of  DD  kernel  energies,  has 
been  shown  to  be  a  transferable  kinetic  energy  functional  which  predicts  a  wide  range  of 
bulk  and  defect  properties  in  Aluminum.  To  address  this  challenge,  we  developed  a  local 
variational  formulation  of  these  extended  kinetic  energy  functionals  through  a  system  of  cou¬ 
pled  Helmholtz  equations  (Radhakrishnan  &  Gavini,  2010).  The  resulting  local  real-space 
saddle-point  formulation  was  discretized  using  a  finite-element  basis,  and  the  computed 
ground-state  properties  were  in  excellent  agreement  with  the  Fourier  space  calculations. 
The  quasi-continuum  reduction  of  orbital-free  DFT  including  the  Dl  and  DD  kinetic  energy 
functionals  was  subsequently  implemented  taking  advantage  of  the  adaptive  coarse-graining 
capability  of  the  finite-element  basis.  We  subsequently  conducted  a  cell-size  study  on  the 
formation  energies  of  vacancies  by  considering  cell-sizes  up  to  a  million  atoms.  These  studies 
have  revealed  larger  cell-size  effects  in  the  computed  energetics  than  previously  reported, 
and  have  demonstrated  that  long-ranged  electronic  and  elastic  fields  play  an  important  role 
in  governing  the  energetics  of  defects  (cf.  figures  1  and  2). 


Figure  1:  Effect  of  cell-size  on  mono¬ 
vacancy  formation  energy.  Fi§ure  2:  Effect  of  cell-size  on  di-vacancy 

binding  energy. 

As  part  of  this  work  we  also  conducted  a  mathematical  analysis  of  the  quasi-continuum  re¬ 
duction  technique  (Gavini  &  Liu,  2011).  In  particular,  by  perturbation  method  and  multiple 
scale  analysis,  we  provided  a  formal  justification  for  the  validity  of  numerical  coarse-graiuiug 
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of  various  fields  in  the  quasi-continuum  reduction  of  orbital-free  DFT  by  taking  the  homog¬ 
enization  limit.  Further,  we  derived  the  homogenized  equations  that  govern  the  behavior  of 
electronic  fields  in  regions  of  smooth  deformations.  Using  Fourier  analysis,  we  determined 
the  far-held  solutions  for  the  electronic  and  clastic  fields  in  the  presence  of  local  defects,  and 
subsequently  estimate  cell-size  effects  in  computed  defect  energies.  Our  studies  have  indi¬ 
cated  that  even  for  weak  defects  like  vacancies,  the  perturbations  in  electronic  fields  in  the 
presence  of  the  defects  have  non-vanishing  contributions  to  the  energy  up  to  cell-size  of  the 
order  of  1000  atoms  which  is  consistent  with  the  results  from  our  numerical  investigations. 


2  Real-space  formulation  of  Kohn-Sham  density  func¬ 
tional  theory  and  Mathematical  analysis: 

The  validity  of  orbital-free  DFT  is  restricted  only  to  material  systems  whose  electronic  struc¬ 
ture  is  close  to  a  free-elcctron  gas — simple  metals  and  Aluminum.  In  order  to  over  come 
this  limitation  we  have  developed  a  local  real-space  variational  formulation  of  Kohn-Sham 
DFT  (Suryanarayana  et  al.,  2010).  In  particular,  the  extended  electrostatic  interactions  are 
reformulated  as  a  local  variational  problem  by  solving  the  Poisson  equation,  and  the  compu¬ 
tation  of  ground-state  properties  in  Kohn-Sham  DFT  reduces  to  a  local  saddle-point  problem 
in  the  vector  containing  the  single-electron  wave-functions  and  the  electrostatic  potential. 
We  note  that  this  variational  formulation  is  equivalent  to  solving  an  eigen-value  problem, 
with  the  eigenfunctions  representing  the  single-electron  wave-functions,  in  a  self-consistent 
iteration  representing  a  fixed  point  iteration  of  the  Kohn-Sham  map.  The  numerical  dis¬ 
cretization  of  the  formulation  has  been  implemented  using  a  finite-element,  basis  constructed 
from  a  finite-element  triangulation  of  the  domain  with  linear  tetrahedral  elements.  Both  the 
non-linear  variational  formulation  and  the  self-consistent  fixed  point  iteration  approach  solv¬ 
ing  the  eigenvalue  problem  have  been  numerically  implemented  and  benchmark  studies  have 
been  conducted.  An  important  outcome  of  the  present  study  was  the  excellent  scalability  of 
the  discrete  formulation  on  parallel  computing  architectures  afforded  by  the  finite-element 
basis  as  well  as  the  ability  to  handle  both  pseudopotential  and  all-electron  calculations  with 
complex  boundary  conditions.  However,  this  study  also  demonstrated  the  computational 
inefficiency  of  linear  finite-element  basis  functions  in  electronic  structure  calculations  per¬ 
formed  using  Kohn-Sham  DFT.  It  was  seen  that  about  100,000  linear  FE  basis  functions 
per  atom  were  required  to  achieve  chemical  accuracy,  which  is  significantly  larger  than  other 
basis  functions  like  plane-wave  basis  and  Gaussian  basis  functions.  In  order  to  address  this 
bottleneck,  we  have  developed  the  framework  for  higher-order  finite-element  discretization 
that  is  discussed  subsequently. 

As  part  of  this  work  we  also  conducted  the  mathematical  analysis  of  the  variational 
formulation  of  Kohn-Sham  DFT  (Suryanarayana  et  ah,  2010).  In  particular,  using  the  di¬ 
rect  method  in  the  calculus  of  variations,  we  proved  the  existence  of  solutions  of  variational 
problem.  A  key  idea  used  in  this  proof  includes  Sobolev  embeddings  to  prove  the  lower  semi¬ 
continuity  in  the  weak  topology.  Further,  we  show  the  T— convergence  of  the  finite-element 
discretization,  which  is  a  variational  notion  of  convergence  for  non-linear  functionals.  As  a 
further  part  of  this  study  we  also  prove  the  T— convergence  of  the  finite-element  discretiza¬ 
tion  with  numerical  quadratures  of  sufficient  accuracy.  A  key  idea  in  this  proof  involves  the 
use  inverse  inequalities  to  show  the  equi-coercivity  of  the  sequence  of  functionals  in  the  weak 
topology. 
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3  Computational  efficiency  afforded  by  higher  order 
finite-element  discretization  in  electronic  structure 
calculations: 

Our  studies  indicated  that  linear  FE  basis  functions  are  computationally  inefficient  in  elec¬ 
tronic  structure  calculations  which  require  chemical  accuracy.  To  this  end  we  explored  the 
computational  efficiency  afforded  by  higher-order  finite-element  discretizations.  In  particu¬ 
lar,  we  developed  a  framework  for  higher  order  finite-element  discretizations  which  include 
quadratic  tetrahedral  FE  basis  functions,  and  hexahedral  node-based  Lagrange  basis  func¬ 
tions  up  to  10th  order  which  includes  spectral  elements.  We  conducted  systematic  studies 
to  assess  the  computational  efficiency  afforded  by  higher-order  FE  discretizations  for  both 
orbital-free  DFT  (Motamarri  et  ah,  2011)  and  Kohn-Sham  DFT  (Motamarri  et  ah,  2012). 
Our  benchmark  numerical  studies  demonstrated  optimal  rates  of  convergence  for  these  dis¬ 
cretization  based  on  linear  finite-element  analysis  although  the  present  problem  is  non-linear. 
Further,  more  interestingly,  the  optimal  rates  of  convergence  were  also  observed  for  all¬ 
electron  calculations,  which  represents  a  system  with  Coulomb  singularity  in  the  electrostatic 
potential.  To  the  best  of  our  knowledge  there  is  no  mathematical  theory  which  guarantees 
the  optimal  rates  of  convergence  in  all-electron  calculations,  and  the  present  studies  provide 
the  first  numerical  evidence  for  optimal  convergence.  While  higher-order  discertizations  are 
more  accurate  and  require  lesser  number  of  basis  functions  to  achieve  chemical  accuracy  in 
comparison  to  a  lower-order  discretization,  the  per  basis  computational  cost  is  also  higher 
for  higher-order  discretizations.  Thus,  in  order  to  assess  the  computational  efficiency  of 
higher  order  finite-element  discretizations  we  conducted  a  study  to  compare  the  CPU  times 
required  to  compute  the  electronic  structure  of  benchmark  problems  using  various  orders 
of  finite-element  basis  functions.  Our  studies  indicate  that  a  4th  order  FE  basis  function 
provides  a  staggering  100-1000  fold  computational  savings  in  comparison  to  a  linear  FE  basis 
function.  This  clearly  demonstrates  the  computational  efficiency  afforded  by  higher  order 
discretizations  in  electronic  structure  calculations.  Further,  we  observed  that  spectral  FE 
basis  functions  improve  the  conditioning  of  the  system  and  result  in  a  faster  convergence  of 
the  iterative  algorithms. 


Figure  3:  Convergence  of  error 
in  ground-state  energy  of  a  35- 
atom  Barium  cluster  for  various 
orders  of  finite-elements. 


Figure  4:  Relative  perfor-  Figure  5:  Electron-density 

mance  of  various  orders  of  finite-  contours  of  an  Aluminum 
elements  measured  in  terms  of  cluster  with  172  atoms, 
computational  time. 


The  spectral  FE  basis  functions  have  a  special  role  in  Kohn-Sham  DFT  calculations.  The 
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discrete  eigenvalue  problem  while  using  spectral  FE  basis  functions  reduces  to  a  standard 
eigenvalue  problem  for  which  efficient  numerical  techniques  have  been  previously  developed. 
In  particular,  in  the  present  work,  we  took  advantage  of  the  Chebyschev  acceleration  tech¬ 
nique  to  solve  for  the  subspace  spanned  by  the  single-electron  wave-functions.  Upon  comput¬ 
ing  the  eigenspace,  we  project  the  Hamiltonian  into  this  eigenspace  to  solve  a  smaller  eigen¬ 
value  problem  in  the  subspace,  which  is  numerically  more  efficient  than  solving  the  eigenvalue 
problem  in  the  original  space.  We  conducted  benchmark  studies  to  compare  the  performance 
of  our  implementation  with  existing  softwares  like  ABINIT  and  Gaussian.  The  proposed  ap¬ 
proach  to  the  solution  of  the  Kohn-Sham  DFT  problem  with  higher-order  FE  basis  functions 
is  comparing  favorably  with  existing  highly  optimized  commercial  codes  (Motamarri  et  ah, 
2012).  Moreover,  we  note  that  the  finite-element  discretization  offers  many  other  advantages, 
which  include  the  ability  to  handle  complex  geometries  and  boundary  conditions,  good  seal- 
ability  on  parallel  computing  architectures,  and  the  ability  to  handle  both  pseudopotential 
and  all-electron  calculations. 

4  Tensor-structured  techniques  for  electronic  structure 
calculations: 


Figure  6:  Electron  density  contours  (computed  ..  , 

.  .  „„„  figure  7:  The  convergence  of  the  approxi- 

trom  Of  Df  1)  ot  a  cluster  containing  odd  Aiu- 

mation  error  with  respect  to  the  Tucker  ten- 
minum  atoms  with  random  replacement  ot  20  . 

..  ,  „  sor  rank  m  the  trobemus  and  energy  norms 

Aluminium  atoms  by  10  Lithium  atoms  and  10 

„  ;  , .  ,  ,  „  tor  varying  univariate  grid-size  n  =  64,  128. 

vacancies.  Contours  at  sectional  plane  y  =  —8. 

In  the  course  of  this  project,  a  new  research  direction  emerged.  Low-rank  tensor  represen¬ 
tations  of  the  multi-dimensional  data  based  on  the  canonical  and  Tucker  models  have  been  in 
use  in  signal  processing  and  computer  science  (De  Lathauwer  et  ah,  2000a, b;  Kolda  &  Bader, 
2007).  Recent  developments  in  tensor-structured  methods  also  dramatically  reduce  the 
complexity  of  the  numerical  evaluation  of  the  multivariate  functions  and  operators  in 
d  >  3.  Algebraically  separable  tensor  approaches  have  recently  been  applied  to  electronic 
structure  calculations  in  the  context  of  Hartree-Fock  equation  (Khoromskij  et  ah,  2011; 
Khoromskaia  et  ah,  2011).  We  conducted  a  study  to  investigate  dependence  of  the  rank 
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of  the  Tucker  approximation  of  the  electronic  structure  of  systems  with  varying  sizes.  To 
this  end,  we  computed  the  electronic  structure  of  varying  size  of  Aluminum  clusters  which 
include  lxlxl  FCC  unit  cells  to  5  x  5x  5  FCC  unit  cells  using  orbital- free  DFT.  The 
Tucker  representation  of  this  data  was  subsequently  computed  and  analyzed  in  collabora¬ 
tion  with  Dr.  Blesgen  and  Dr.  Khoromskaia  at  the  Max-Planck  Institute  for  Mathematical 
Sciences.  The  analysis  showed  that  a  rank  of  less  than  10  was  sufficient  to  approximate  the 
electronic  structure  of  the  various  systems  with  varying  sizes,  and  the  rank  was  independent 
of  the  number  of  atoms  in  the  system  for  the  systems  considered  (Blesgen  et  ah,  2011).  As  a 
more  stringent  test,  we  replaced  at  random  20  Aluminum  atoms  in  a  non-periodic  5x5x5 
FCC  Aluminum  cluster  with  666atms  with  with  10  vacancies  and  10  Lithium  atoms.  Our 
studies  show  that  a  Tucker  rank  of  20  was  sufficient  to  approximate  the  electron-density  of 
this  materials  system,  suggesting  that  tensor-structured  techniques  based  on  Tucker  repre¬ 
sentation  may  be  efficient  representations  of  electronic  structures.  These  studies  open  the 
possibility  of  computing  both  the  orbital-free  DFT  problem  and  Kohn-Sham  DFT  problem 
using  tensor-structured  techniques,  and  the  development  of  these  algorithms  is  currently 
being  pursued. 

5  Application  studies: 


(a)  (b) 

Figure  8:  (a)  Dependence  of  vacancy  formation  energy  on  macroscopic  deformations  com¬ 
puted  using  quasi-continuum  orbital-free  DFT.  (b)  Contours  of  perturbation  in  electron- 
density  in  the  presence  of  a  vacancy;  left:  no  imposed  deformation;  right:  0.33  volumetric 
strain. 


We  used  the  quasi-continuum  reduction  of  orbital-free  DFT  to  investigate  if  the  defect- 
core  governed  by  quantum-mechanical  interactions  plays  a  crucial  role  in  the  energetics  of 
defects  by  studying  vacancies  in  Aluminum.  We  find  that  vacancy  core-energy  is  signif¬ 
icantly  influenced  by  the  state  of  deformation  at  the  vacancy-core,  especially  volumetric 
strains  (Gavini,  2009).  Insights  from  the  core  electronic  structure  and  computed  displace¬ 
ment  fields  show  that  this  dependence  on  volumetric  strains  is  closely  related  to  the  changing 
nature  of  the  core-structure  under  volumetric  deformations  (cf.  figure  8).  These  results  are 
in  sharp  contrast  to  mechanics  descriptions  based  on  elastic  interactions  that  often  consider 
defect  core-energies  as  an  inconsequential  constant.  Calculations  suggest  that  the  variation 
in  core-energies  with  changing  macroscopic  deformations  is  quantitatively  more  significant 
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than  the  corresponding  variation  in  relaxation  energies  associated  with  clastic  fields.  Upon 
studying  the  influence  of  various  macroscopic  deformations,  which  include  volumetric,  uni¬ 
axial,  biaxial  and  shear  deformations,  on  the  formation  energies  of  vacancies,  we  show  that 
volumetric  deformations  play  a  dominant  role  in  governing  the  energetics  of  these  defects. 
Further,  by  plotting  formation  energies  of  vacancies  and  di- vacancies  against  the  volumetric 
strain  corresponding  to  any  macroscopic  deformation,  we  find  that  all  variations  in  the  for¬ 
mation  energies  collapse  on  to  a  universal  curve.  This  suggests  a  universal  role  of  volumetric 
strains  in  the  energetics  of  vacancies.  Implications  of  these  results  in  the  context  of  dynamic 
failure  in  metals  due  to  spalling  are  analyzed. 
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